package scu.maqiang.eigenvalue;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES1L2PE;
import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.EigenSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class PEEigenValue1DL21 {

	public static void main(String[] args) {
		Mesh1L2 mesh = new Mesh1L2();
		mesh.line(0, 1.0, 4096);
		mesh.translate(1.0);
		//mesh.scale(Math.PI);
		
		FES1L2PE fs = new FES1L2PE(mesh);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		SRMatrix M = new SRMatrix(fs.GetNdof());
		
		int N = 10;
		double[] D = new double[N];
		double[][] V = new double[N][fs.GetNdof()];
		double[] x = new double[fs.GetNdof()];
		Arrays.fill(x, 1.0);
		
//		double G1 = 1.0, ep1 = 0.5, e1 = 0.1;
//		fs.assemblePEStiff(new double[] {G1, ep1, e1}, BVPType.COMMON, A);
//		fs.assemblePEMass(new double[] {1.0}, BVPType.CONSIST_MATRIX, M);
		
		double G1 = 1.0, G2 = 2.0, lambda = 0.5, ep = 0.2, e1 = 0.4, e2 = 0.3;
		fs.assemblePEStiff(new double[] {G1, G2, lambda, ep, e1, e2}, BVPType.AXISYMMETRIC, A);
		fs.assemblePEMass(new double[] {1.0}, BVPType.AXIS_CONSIST_MASS, M);
		//M.sort();
		//System.out.println(M);
		fs.applyBC_MBN(A, x, Direct.X, 0.0, 1);
		fs.applyBC_MBN(A, x, Direct.Y, 0.0, 2);
		//fs.applyBCZero(A, x, Direct.X, BCType.RCE, 1);
		//fs.applyBCZero(A, x, Direct.Y, BCType.RCE, 2);
		//fs.applyBCZero(A, x, Direct.X, BCType.RCE, 1, 2);
		//System.out.println(A);
		EigenSSolver solver = new EigenSSolver();
		solver.inversePowerTran_PLU(A, M, fs.getnBand(), x, 20.0, N, D, V);
		System.out.println(MVO.toString(D));
		//mesh.toTecplot("Eigenvalues1D.dat", V);
		String outputFilePreFix = "AxisPEEigenValue";
		double[][] u = new double[2][mesh.getNv()];
		for(int i = 0; i < N; i++) {
			fs.extractComponent(V[i], u);
			mesh.toTecplot(outputFilePreFix + (i+1) + ".dat", u);
		}
	}

}
